#' Derive reflectance and temperature from Landsat
#'
#' Landsat ETM+ raw data is formatted as a digitial number. Several calculations must be performed to convert the digital number into temperature and reflectance.
#'
#' @name LReTe-package
#' @aliases LReTe
#' @docType package
#' @title Landsat ETM+ reflectance and temperature functions.
#' @author Joseph Henry \email{joseph.henry@@sydney.edu.au} and Willem Vervoort \email{willem.vervoort@@sydney.edu.au}
#' @references NASA Handbook (to add)
#' @keywords package, Reflectance, Radiance

### CONDITION OF USE ###
# Joseph reserves the right to use this code or its derivatives in the creation of an R package.

# Temporary constants
	DATE.FORMAT.3 <- "%Y%m%d"
	DATE.FORMAT.5 <- "%Y-%m-%d"
	DATE.FORMAT.6 <- "%j"
	SYMBOL.EQUAL <- "="
	SPACE.SINGLE <- " "
	SPACE.NONE <- ""
	
							
	shad_cell_func <- function(input.a,input.b){
		shadow_cells <- cellFromRowCol(input.a,input.b[,1],input.b[,2])
		return(shadow_cells)
	
	}
	L7CLIP <- function (x, y, bbox){
		#' @param x input
		#' @param y output
		#' @param bbox (ul, ur, ll, lr (new extent))
		r <- raster (x)
		e <- extent(bbox)
		rc <- crop(r, e)
		rm (r); rm (e)
		writeRaster(rc, y)
		rm (rc)
	}
	L7ConstantBit <-  function (x,y) {
		output <- substring (x, y, nchar (x))
		return (output)
	}
	QC_information_func <-  function (x){
		Numbers_and_constants <-  gsub (" ", "", unlist (strsplit (x, "="))) # ?	
		R1.SELECT <-  which (nchar (Numbers_and_constants) <= 5) # ?
		R1 <-  Numbers_and_constants[R1.SELECT]
		R2 <-  matrix (R1, ncol = 2, byrow = TRUE) # Convert to matrx
		return (R2) # Return value
		rm (R1);rm (R2); rm (Numbers_and_constants); rm (R1.SELECT) # Clear memory 
	}
	L7GainFromFN <-  function (x, y){
		output <- substring (x, nchar (x) - y, nchar (x))
	}
	L7GetAcquisitionDate  <-  function (x, y) {
		#' Get the tile date from the file name
		#' @param x Metadata filename (character) 
		#' @param y number characters from the end (numeric)
		#' @return date (character)	
		output <-	gsub (" ", "", substring (x, nchar (x) - y, nchar (x)))
	}
	L7.MetadataSelectRow <- function(x){
		#' Select the appropriate row in the metadata file. File lengths differ over span of the data.
		#' @param x Scanned Landsat tile metadata
		#' @return row list 
		# Predefined data
		 L7.META.LENGTH.OPTIONS <- c (173, 175, 176) # Possible metadata lengths
		 L7.META.ACUISITIONDATE.ROW.OPTIONS <- list (20,21) # Options for row containing acuisition date data
		 L7.META.QC.ROW.OPTIONS <-  list (81:98, 83:100) # Options for rows containing quality control data
		 L7.META.GAIN.ROW.OPTIONS <-  list (111:119, 113:121) # Options for rows containing gain data
		 L7.META.SUN.ROW.OPTIONS <-  list (138:139, 140:141) # Options for rows containing sun angle data
		# Extract information about input
		 x_length <- length(x) # Metadata length
		# Length test
		 if(x_length == L7.META.LENGTH.OPTIONS[1]){
				L7.META.ACUISITIONDATE.ROW <- L7.META.ACUISITIONDATE.ROW.OPTIONS[[1]]
				L7.META.QC.ROW <- L7.META.QC.ROW.OPTIONS[[1]]
				L7.META.GAIN.ROW <- L7.META.GAIN.ROW.OPTIONS[[1]]
				L7.META.SUN.ROW <- L7.META.SUN.ROW.OPTIONS[[1]]	
			}
			if(x_length == L7.META.LENGTH.OPTIONS[2] | x_length ==L7.META.LENGTH.OPTIONS[3]){
				L7.META.ACUISITIONDATE.ROW <- L7.META.ACUISITIONDATE.ROW.OPTIONS[[2]]
				L7.META.QC.ROW <- L7.META.QC.ROW.OPTIONS[[2]]
				L7.META.GAIN.ROW <- L7.META.GAIN.ROW.OPTIONS[[2]]
				L7.META.SUN.ROW <- L7.META.SUN.ROW.OPTIONS[[2]]	
			}
		# Length test output
		 MetadataSelectRow.LIST.OUT <- list (L7.META.ACUISITIONDATE.ROW,L7.META.QC.ROW, L7.META.GAIN.ROW, L7.META.SUN.ROW) # List results for output
	 return (MetadataSelectRow.LIST.OUT) # Return
	}
	L7SunInformation <-  function (MTL_all, y){
		#' Extract sun information from the metadata
		#' @param MTL_all scanned metadata file
		#' @param y return option where azimuth and elevation are 1 or 2, respectfully (numeric)
		#' @return Sun information (numeric) 
		# Determine metadata length characteristics
		 MTL_ALL_CHARACTERISTICS <- L7.MetadataSelectRow (MTL_all)  
		 L7.META.SUN.ROW <- MTL_ALL_CHARACTERISTICS[[4]]
		# Extract relevant metadata
		 RELEVANT_METADATA <- MTL_all[L7.META.SUN.ROW]	
		 Numbers_and_constants <-  (gsub (SPACE.SINGLE, SPACE.NONE, unlist (strsplit (RELEVANT_METADATA, SYMBOL.EQUAL))))
		 R1 <- matrix (Numbers_and_constants, ncol = 2, byrow = TRUE)
		# Output 
		 R2 <- as.numeric (R1 [y,2] )
		 return (R2)
		# Clear memory 
		 rm (R1); rm (R2); rm (Numbers_and_constants)
	}
	FilenamePaste <- function (x, y) {
		#' Filename paste function
		#' @param x names of each Landsat band
		#' @param y path to folder
		#' @return full path of each Landsat band 
		Output <-	paste (y ,paste (x, ".tif", sep = ""), sep = "/")
		# y =  DATA.DIR,"processed/gis/raster/2d/L7_derived/original_projection/extent"
		return (Output)
	}	
	LReTe <- function (L7.BANDNAMES.12345617, OUTPUT.NAME, MTL_all, OR.CHOICE) {
		#' Produce Landsat reflectance and temperature rasters.
		#' @param MTL_all scanned metadata file
		#' @param L7.BANDNAMES.12345617 file path of each Landsat input band
		#' @param OUTPUT.NAME file path for each Landsat band
		#' @param L7.LIST.SUN
		#' @param OR.CHOICE choice to overwrite rasters (boolean)
		#' @return NULL  (rasters with radiance and temperature in degrees Celcius)
		
		 ## Load preparation
			 L7.PASS.CHOICE <- L7.BANDNAMES.12345617 # Bands for each Landsat tile. Each row is a different tile.
			# Pass raster filenames
			 b1.dir <-  (L7.PASS.CHOICE[1]) # Band 1 path
			 b2.dir <-  (L7.PASS.CHOICE[2]) # Band 2 path
			 L7.DIR.PASSCHOICE.BAND3 <-  (L7.PASS.CHOICE[3]) # Band 3 path
			 L7.DIR.PASSCHOICE.BAND4 <-  (L7.PASS.CHOICE[4]) # Band 4 path
			 b5.dir <-  (L7.PASS.CHOICE[5]) # Band 5 path	
			 b6.dir <-  (L7.PASS.CHOICE[6]) # Band 6 path
			 b7.dir <-  (L7.PASS.CHOICE[7]) # Band 7 path
			# Assign Band names	
			 Band.name <- paste ("Band_", seq (1, 8, 1), sep = "")
			 Filter.set.1 <- paste ("Filter_p1_", seq (1, 8, 1), sep = "")

			## Load data		
			 L7.RESULTS.METADATAROW <- L7.MetadataSelectRow(MTL_all) # Get correct rows
			 # Assign values
				 L7.META.ACUISITIONDATE.ROW <- L7.RESULTS.METADATAROW[[1]]
				 L7.META.QC.ROW <- L7.RESULTS.METADATAROW[[2]]
				 L7.META.GAIN.ROW <- L7.RESULTS.METADATAROW[[3]]
				 L7.META.SUN.ROW <- L7.RESULTS.METADATAROW[[4]]
					# Satellite gain information
						L7.GAININFORMATION <- gsub (" ", "", L7GainFromFN (gsub ("  ", "", MTL_all[L7.META.GAIN.ROW]), 1))	
					# Image acuisition date					
						L7.ACUISITIONDATE <- L7GetAcquisitionDate (MTL_all[L7.META.ACUISITIONDATE.ROW], 10)
					# Spectral radiance choice dependent upon acquisition date
						if(as.Date (L7.ACUISITIONDATE, DATE.FORMAT.5) > as.Date ("20000701", DATE.FORMAT.3)){
								L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE.RAW <- L7.READ.LUT.SPECTRALRADIANCERANGE[3 : nrow (L7.READ.LUT.SPECTRALRADIANCERANGE), c (1, 6 : 9)]
								L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE <- rbind(L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE.RAW[1:5,],L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE.RAW[6,],L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE.RAW[6,],			L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE.RAW[7:8,])	# 6 Repeated because there are two band 6 values
						}		
					BAND <- L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE[,1]
					L7.QCAL <- QC_information_func(MTL_all[L7.META.QC.ROW])
					L7.QCAL.MAX <- L7.QCAL[,1]
					L7.QCAL.MIN <- L7.QCAL[,2]
					LMIN <- c(); LMAX <- c()	# Establish vector
				 # Choice for if use high or low gain
					for(i in seq(L7.GAININFORMATION)){
					 if (L7.GAININFORMATION [i] == "H"){
					 LMIN[i] <- L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE[i,4]
					 LMAX[i] <- L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE[i,5]
					 }
					 if (L7.GAININFORMATION[i] == "L"){
					 LMIN[i] <- L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE[i,2]
					 LMAX[i] <- L7.READ.LUT.SPECTRALRADIANCERANGE.CHOICE[i,3]
					 }
					}
				 # Metadata for Landsat tile		
				  L7.LIST.METADATA <- list (data.frame (BAND, 
					L7.WORD.PART.1 = (as.numeric (LMAX) - as.numeric (LMIN)) / (as.numeric (L7.QCAL.MAX) - as.numeric(L7.QCAL.MIN)), 
					L7.WORD.QCAL.MIN = as.numeric(L7.QCAL.MIN), L7.WORD.LMIN = as.numeric(LMIN))
					)
			 names(L7.LIST.METADATA) <- L7.FN.PASSNAME.CONSTANTBITEXTRACT			
				L7.RADIANCE.CHOICE <- data.frame (L7.LIST.METADATA[[1]])
				L7.SUN.ELEVATION <- L7SunInformation (MTL_all, 2)			
				L7.SUN.AZIMUTH <- L7SunInformation (MTL_all, 1)
			 # Rasters
				Band_1 <- raster (b1.dir) # Out by itself here because it is used as a template
			
			## Clean data	
			 # Sun azimuth/elevation information
				L7.SUN.ELEVATION.TAN.DEGREES <- 90 - L7.SUN.ELEVATION # Solar zenith angle (reciprocal of sun elevation angle)
				L7.SUN.ELEVATION.TAN.RADIANS <- L7.SUN.ELEVATION.TAN.DEGREES * ( pi / 180)
				L7.SUN.AZIMUTH.TAN.DEGREES <- L7.SUN.AZIMUTH - 90
				L7.SUN.AZIMUTH.TAN.RADIANS <- L7.SUN.AZIMUTH.TAN.DEGREES * ( pi / 180)

			## Do: Convert Bands 1, 2, 3, 4, 5 and 6 to reflectance
			 L7.COEFFICIENT.ZHU <- L7.COEFFICIENT.ZHU # Multiplyer used by Zhu. Maybe this is used to convert from floating poin to integer?
			 L7.LOOP.BANDS34 <- 1:7 # Loop for each band
			 # Loop through and assign values for each Band	
				for (i in 1 : ( length(L7.LOOP.BANDS34))){
				 gc();
				 print (paste(i,"---",Sys.time ()));
				 if (i <= 6){ # Band 6 comes in two forms so one file has to skipped. Therefore, Band 7 equals 8, not 7.
					ref.row <- L7.LOOP.BANDS34 [i]
				 }
				 if (i == 7){
					ref.row <- 8
				 }
				 # Vectorise				
					if(L7.LOOP.BANDS34[i] == 1){					
					 Band_1_vec_raw <- Band_1[] # Convert rasters to vectors
					 L7.REFLECTANCEVALUE <- Band_1_vec_raw
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE # Convert NA back to L7.FILLVALUE
					 rm (Band_1_vec_raw)
					}
					if(L7.LOOP.BANDS34[i] == 2){					
					 Band_2 <- raster (b2.dir)
					 Band_2_vec_raw <- Band_2[]
					 rm (Band_2);
					 L7.REFLECTANCEVALUE <- Band_2_vec_raw
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (Band_2_vec_raw)
					}				
					if(L7.LOOP.BANDS34[i] == 3){		
					 L7.RASTER.BAND.3 <- raster (L7.DIR.PASSCHOICE.BAND3)
					 L7.VECTOR.BAND.3.RAW <- L7.RASTER.BAND.3[]
					 rm (L7.RASTER.BAND.3); 	
					 L7.REFLECTANCEVALUE <- L7.VECTOR.BAND.3.RAW
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (L7.VECTOR.BAND.3.RAW)
					}
					if(L7.LOOP.BANDS34[i] == 4){	
					 L7.RASTER.BAND.4 <- raster (L7.DIR.PASSCHOICE.BAND4)
					 L7.VECTOR.BAND.4.RAW <- L7.RASTER.BAND.4[]
					 rm (L7.RASTER.BAND.4);					
					 L7.REFLECTANCEVALUE <- L7.VECTOR.BAND.4.RAW
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (L7.VECTOR.BAND.4.RAW)
					}
					if(L7.LOOP.BANDS34[i]==5){		
					 Band_5 <- raster (b5.dir)
					 Band_5_vec_raw <- Band_5[]
					 rm (Band_5);
					 L7.REFLECTANCEVALUE <- Band_5_vec_raw
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (Band_5_vec_raw)
					}
					if(L7.LOOP.BANDS34[i]==6){		
					 Band_6 <- raster (b6.dir) # Low gain Band 6 used. However, the other one can be used. I don't understand the impact of using the other one (JH)
					 Band_6_vec_raw <- Band_6[]
					 rm (Band_6); 					
					 L7.REFLECTANCEVALUE <- Band_6_vec_raw
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (Band_6_vec_raw)
					}
					if(L7.LOOP.BANDS34[i]==7){			
					 Band_7 <- raster (b7.dir)
					 Band_7_vec_raw <- Band_7[]
					 rm (Band_7);					
					 L7.REFLECTANCEVALUE <- Band_7_vec_raw
					 gc();
					 L7.REFLECTANCEVALUE[is.na (L7.REFLECTANCEVALUE) == TRUE] <- L7.FILLVALUE
					 rm (Band_7_vec_raw)
					}
					# Derive radiance
					 L7.RADIANCEVALUE <- (L7.RADIANCE.CHOICE[ref.row,2]) * (L7.REFLECTANCEVALUE - L7.RADIANCE.CHOICE[L7.LOOP.BANDS34[i],3]) +L7.RADIANCE.CHOICE[L7.LOOP.BANDS34[i],4]
					 rm (L7.REFLECTANCEVALUE)
					# Derive reflectance 
					 if(L7.LOOP.BANDS34[i] != 6){ # Band 6 is for temperature.
					  L7.ACQUISITION.JULIANDAY <- as.numeric (format (strptime (L7.ACUISITIONDATE,DATE.FORMAT.5), DATE.FORMAT.6))
					doy_app <- which(L7.READ.LUT.EARTHSUNDISTANCE[,1] %in% L7.ACQUISITION.JULIANDAY)
					ES_DISTANCE <- L7.READ.LUT.EARTHSUNDISTANCE[doy_app,2]
						# Cosine radians used because will only produce positive reflectance, unlike cosine degrees
						Reflectance <- (pi * L7.RADIANCEVALUE * ES_DISTANCE^2) / (L7.READ.LUT.SOLARSPECTRAIRRADIANCE[L7.READ.LUT.SOLARSPECTRAIRRADIANCE[,1] == L7.RADIANCE.CHOICE[ref.row,1],2] * cos (L7.SUN.AZIMUTH.TAN.RADIANS))
						rm(L7.RADIANCEVALUE)
					 }
					 
					# Output data
					 if(L7.LOOP.BANDS34[i] == 1){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						gc();
						rm (Band_vec)
						writeRaster(BAND.OUT, OUTPUT.NAME[1],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 2){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						gc();
						rm (Band_vec)
						writeRaster(BAND.OUT, OUTPUT.NAME[2],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 3){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						gc();
						rm (Band_vec)
						writeRaster(BAND.OUT, OUTPUT.NAME[3],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 4){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						gc();
						rm (Band_vec)
						writeRaster(BAND.OUT, OUTPUT.NAME[4],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 5){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						gc();
						rm (Band_vec)
						writeRaster(BAND.OUT, OUTPUT.NAME[5],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 6){
						gc();
						Band_6_vec_rad <- L7.RADIANCEVALUE
						rm (L7.RADIANCEVALUE)
						## Convert Band 6 to at sensor temperature
						K2 <- 1282.71 # Kelvin calibration constant
						K1 <- 666.09 # Watts...calibration constant
						gc ()
						Band_6_vec <- K2 / log ((K1 / Band_6_vec_rad) + 1)
						gc ()
						Band_6_vec_Cel <- 100 * (Band_6_vec - 273.15) # Convert Band 6 Kelvin to Celcius with 0.01 scale factor
						rm (Band_6_vec)
						gc()
						writeRaster(setValues(Band_1, Band_6_vec_Cel), OUTPUT.NAME[6],overwrite = OR.CHOICE) # Output
						rm (Band_6_vec_Cel) # Remove data
					 }
					 if(L7.LOOP.BANDS34[i] == 7){
						gc();
						Band_vec <- Reflectance * L7.COEFFICIENT.ZHU
						rm (Reflectance)
						BAND.OUT <- setValues(Band_1, Band_vec)
						rm (Band_vec)
						gc();
						writeRaster(BAND.OUT, OUTPUT.NAME[7],overwrite = OR.CHOICE) # Output
						rm (BAND.OUT) # Remove data
					 }
				}
	# Remove data
		rm (Band_1);
		return (c (L7.SUN.AZIMUTH.TAN.RADIANS,L7.SUN.AZIMUTH,L7.SUN.ELEVATION.TAN.RADIANS ))
	}			


### Nice function to check memory
# NOTE: This is copied from a stackoverflow thread.
# Source: http://stackoverflow.com/questions/1358003/tricks-to-manage-the-available-memory-in-an-r-session
.ls.objects <- function (pos = 1, pattern, order.by,
					decreasing=FALSE, head=FALSE, n=5) {
napply <- function(names, fn) sapply(names, function(x)
									 fn(get(x, pos = pos)))
names <- ls(pos = pos, pattern = pattern)
obj.class <- napply(names, function(x) as.character(class(x))[1])
obj.mode <- napply(names, mode)
obj.type <- ifelse(is.na(obj.class), obj.mode, obj.class)
obj.size <- napply(names, object.size)
obj.prettysize <- sapply(obj.size, function(r) prettyNum(r, big.mark = ",") )
obj.dim <- t(napply(names, function(x)
					as.numeric(dim(x))[1:2]))
vec <- is.na(obj.dim)[, 1] & (obj.type != "function")
obj.dim[vec, 1] <- napply(names, length)[vec]
out <- data.frame(obj.type, obj.size,obj.prettysize, obj.dim)
names(out) <- c("Type", "Size", "PrettySize", "Rows", "Columns")
if (!missing(order.by))
	out <- out[order(out[[order.by]], decreasing=decreasing), ]
	out <- out[c("Type", "PrettySize", "Rows", "Columns")]
	names(out) <- c("Type", "Size", "Rows", "Columns")
if (head)
	out <- head(out, n)
out
}
